function [par_long_opt,value_opt,aic,prob_opt,distrib_opt,EQD2_opt,Effect_opt,vector_loglike] = MaxLikelihood_TCP(data_TCPexp,par_long0,ind_opt0, T0, dT, n_loop,logistic,dep_a,dep_b,treatment_type)  
%  This code performs fits to RT dose response data in terms of tumor control
%  probability (TCP) using clinical information from a database: dose per 
%  fraction, number of RT fractions, treatment time, and tumor control 
%  probability and number of patients/metastases treated. Two different TCP
%  formulations (Logit or Poisson) and several expressions of the surviving 
%  fraction (SF) models can be used with functional dependencies (dep_a, dep_b):

% log (SF)=n.*d.*alpha.*(1+a.*d.^dep_a+d.*(1+b.*d.^dep_b)/alphabeta)-lambda*max(0,t_trat-tk)
% The LQL model is run when dep_a=-1, with: 
% log(SF)=n.*alpha.*(d+2.*(a.*d+exp(-a.*d)-1)./(alphabeta.*a^2))-lambda*max(0,t_trat-tk)

% A simulated annealing optimization algorithm is used to find the
% parameter values that maximize the loglikelihood associated to the fits.

%Input:
%        1. data_TCPexp: database array with clinical RT treatment and treatment outcome information 
%            from several patient cohorts. 
%        2. par_long0: TCP model initial parameter values
%              For Poisson TCP formulation: par_long0=[alpha alphabeta a b g lambda_prolif t_kickoff TCPmax]
%              For Logit TCP formulation: par_long0=[alpha/beta a b D50 gamma lambda_prolif t_kickoff TCPmax] 
%        3: ind_opt: indices of the par_long0 vector corresponding to parameters free to change during the optimization process.
%        4. T0: Initial temperature of the simulated annealing optimization algorithm
%        5. dT: Cooling rate of the simulated annealing optimization algorithm
%        6. n_loop: number of iterations performed in the simulated annealing optimization algorithm
%        7. logistic: binary input specifying the TCP formulation (1 for logistic, 0 for Poissonian)  
%        8. dep_a: functional dependence with dose of the alpha term in the expression of the surviving fraction
%        9. dep_b: functional dependence with dose of the beta term in the expression of the surviving fraction
%        10. treatment_type: String with treatment site ('brain' for brain metastases and 'lung_' for NSCLC.


%Output: 
%          1. par_long_opt: best fitting model parameter values 
%          2. value opt: loglikelihood of the best fit
%          3. aic: best fits AICc values (AICc=Akaike information criterion with sample size correction) 
%          4. prob_opt: best fits probability (use of binomial distribution) 
%          5. distrib_opt: Best fitting TCP model values associated to each cohort of the database
%          6. EQD2_opt: EQD associated to each cohort of the database, calculated with each model functional dependency and parameter value
%          7. Effect_opt: Effect,log(SF), associated to each cohort of the database, calculated with each model functional dependency and parameter value
%          8. vector_loglike: vector with the loglikelihoods associated to the fits during the simulated annealing algorithm convergence


N=1e5; %Number of initial tumor cells (only used in the Poissonian TCP formulation)

%Take treatment information from the database array:
schedules=data_TCPexp(1,7); %up to two adjuvant treatment schedules can be considered (some brain treatment cohorts treated were in this way)
d=data_TCPexp(:,7+2:2:7+2*schedules); % average dose per fraction
inds=find(data_TCPexp(:,1)>170); %Cohorts with more than 170 pts are registered to perform alternative calculation of fitting probability, avoidin divergence of a factorial.
vector_loglike=cell(length(dep_a),1);

%Perform k model optimizations with k functional dependencies specified in dep_a for the alpha term and dep_b for the beta term.
  for k=1:length(dep_a) 
     ind_opt=ind_opt0;
     par_long=par_long0;
     
%Initial a and b parameter value definition depending on the SF functional depencency (lines 44 to 68)
%Thresholds for a and b parameter values are established to avoid negative contributions to treatment effect in the cohort database
     threshold_a=-1/(max(d(:)))^dep_a(k); 
     threshold_b=-1/(max(d(:)))^dep_b(k);
     epsilon=abs(threshold_b/2);
        if logistic==1
            inda=2;
            indb=3; 
        else
            inda=3;
            indb=4;
        end
        
        if dep_a(k)==-1 %LQL condition: a is positive an keeps its sign, therefore a> threshold always
               threshold_a=1e-5; %set to this value to avoid numerical oscilations of LQL function
               par_long(inda)=0.1; %Experimental values between 0.02 and 0.2 (Guerrero et al, 2004 Phys. Med. Biol. 49 4825)         
            elseif dep_a(k)~=0 %Initial condition for parameter a of modified LQ model is set using its threshold value
               par_long(inda)=threshold_a/2;   
            else % If dep_a=0, there is no dependence with a, set a=0 and delete inda from the optimized parameters index list
                par_long(inda)=0;
                ind_opt(ind_opt==inda)=[]; 
        end
        
        if dep_b(k)==0 %If dep_p=0 there is no dependence with b, set b=0 and delete indb from the optimized parameters index list
                par_long(indb)=0;
                ind_opt(ind_opt==indb)=[]; 
        end
%Calculate initial TCP model and loglikelihood associated to the fit with the initial parameter values
            if logistic==1 
                [value,TCPmodel,EQD2_opt(k,:),Effect_opt(k,:),prob_opt(k,:)] = loglikelihood_logistic(par_long,data_TCPexp,inds,dep_a(k),dep_b(k));
            else
                [value,TCPmodel,EQD2_opt(k,:),Effect_opt(k,:),prob_opt(k,:)] = loglikelihood_poisson(par_long,data_TCPexp,inds,N,dep_a(k),dep_b(k),treatment_type);
            end
     value(isinf(value))=1e300; %if prob is too low, the -Log(prod(prob) is too large (a NaN)), to allow optimization, we set it to a large value
     value_opt(k)=value;
     distrib_opt(k,:)=TCPmodel;
     vector_loglike_int=[]; 
     par_opt=par_long(ind_opt);
     par=par_opt;
     ind=length(ind_opt);
     T=T0;
     counter=0;
%Simulated annealing optimization algorithm (lines 92 to 134)
     for i=1:n_loop  
            for j=1:100 
                counter=counter+1;
                %Parameter perturbation (lines 97 to 107)
                %Choose the index of one of the free parameters
                ind2=floor(rand*ind)+1; 
                par2=par;
                  if (ind_opt(ind2)==inda) % Parameter a perturbations keep the sign of the initial value 
                      par2(ind2)=max(threshold_a,par(ind2)+par(ind2)*((rand-0.5)*T/T0));              
                  elseif (ind_opt(ind2)==indb) %The value of parameter b is allowed to change sign
                      par2(ind2)=max(threshold_b,par(ind2)+(T/T0)^1.5*epsilon*(rand-0.5));                    
                  else % The remaining parameter values keep their initial sign
                      %parameter value perturbations of +/- 25% *T/T0  
                      par2(ind2)=par(ind2)+par(ind2)*((rand-0.5)*T/T0);
                   end                
                par_long(ind_opt)=par2; 
                
                %Calculate TCP model (distrib), fit loglikelihood (value2) and other variables (EQD2, Effect and prob) associated to the fit with the perturbed parameter values
                  if logistic==1
                     [value2,distrib,EQD2,Effect,prob] = loglikelihood_logistic(par_long,data_TCPexp,inds,dep_a(k),dep_b(k));
                  else                    
                     [value2,distrib,EQD2,Effect,prob] = loglikelihood_poisson(par_long,data_TCPexp,inds,N,dep_a(k),dep_b(k),treatment_type);
                  end                
                %Oracle function allows survival of nonoptimum param values to avoid trapping 
                surv=oracle(value2-value, T);                   
                  if surv==1               
                      par=par2;                          
                      value=value2;  
                                      
                  end                                
                  %Save parameter value as optimum if AIC is lower than the previously calculated best one      
                  if value<value_opt(k)                 
                      value_opt(k)=value;               
                      par_opt=par;             
                      distrib_opt(k,:)=distrib;             
                      prob_opt(k,:)=prob;   
                      vector_loglike_int(counter)=value;  
                      EQD2_opt(k,:)=EQD2;                                   
                      Effect_opt(k,:)=Effect; 
                 end                               
            end
       T=T*dT;
     end
     par_long(ind_opt)=par_opt;
     par_long_opt(k,:)=par_long;
     vector_loglike{k}=vector_loglike_int;
     length(ind_opt);
     %calculate AICc
     aic(k)=2*value_opt(k)+2*length(ind_opt)+2*length(ind_opt)*(length(ind_opt)+1)/(length(data_TCPexp)-length(ind_opt)-1);
  end
end



function num = oracle(d1,T)  
%This function tests the change in the value of the objective function
%(d1=value2-value1) to check whether a better solution has been found 
%Survival of non optimum solutions is allowed to avoid trapping. The 
%solution survival probability (prob) is proportional to T and inversely 
%proportional to the increase of the objective function (worsening of the solution). 
%T and therefore prob decrease with the number of iterations.

     if d1<0
         num=1;
     else
         
         prob=exp(-d1/T);
         a=rand;
         if a<prob
             num=1;
         else
             num=0;
         end
     end 
end



